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ABSTRACT 

A comparison of several techniques is presented for 
determining upper confidence levels for a system failure 
rate. A series system of components with exponential 
failure rates is examined. Classical computational tech- 
niques are compared with Bayesian techniques in determining 
the upper confidence level of a system failure rate. A 
sensitivity analysis is conducted on several of the 
parameters as part of the comparison. 
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I. INTRODUCTION 



A. BACKGROUND 

Numerous "classic” techniques have been used to compute 
estimates of failure rate and mean time of failure. From 
these estimates standard accepted procedures can be applied 
to establish upper confidence levels (UCL) for the failure 
rate or lower confidence levels (LCL) for the reliability. 

One well established "classic" procedure is to utilize 
the computational methods set forth in Ref. 3 to determine 
upper confidence levels on system failure rates. This is 
the procedure used for all "Classic" and "Semi-Classic" 
methods presented in this paper. 

The application of a "Bayesian" technique may be intui- 
tively appealing to some. Results from previous experiments 
and testing could be applied apriori to current testing 
programs to determine failure rate and reliability. By 
using a prior in such a manner perhaps total system test 
time , number of component failures, etc., could be reduced, 
thereby reducing the overall expense involved with system 
life testing. This would appear to be most appealing when 
testing systems that are extremely expensive. 

The intent of this thesis is to attempt to determine if 
it would be more advantageous to use some "Bayesian" tech- 
nique rather than a more traditional "classic" technique 
to compute upper confidence levels on a system failure rate 
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when the "prior" for the Bayesian method could be chosen 
as "optimistically" as one would desire. 

B. SYSTEM ASSUMPTIONS 

Since it is the technique of computation of the upper 
confidence level that is being investigated, the system 
that is being modeled can be kept simple, yet still 
realistic . 

Even the most complicated system can be broken down to 
a system of components, connected in series, which the 
total failure of any one component will result in the mis- 
sion failure of the system. 

Each type component in this series system experiences 
exponential failure rate. The failure of each type of 
component is assumed to be independent of the failure of 
any other type of component. Wearin and wearout are neg- 
lected and the failure rate is assumed to be constant. 
Components are assumed not to ever be "stillborn." To 
keep the system as simple as possible, each type component 

exhibits an identical failure rate. A, = A.: i = 2 . . . k, 

’ll’ ’ 

k equaling the number of type components that are connected 
in the series system. The following formulas concerning 
system failure rate and system reliability are assumed to 
be valid for this system: 









I m- 



I 



• » 








or R^ = exp 

In keeping the system as simple as possible the 
optimistic "prior" for the Bayesian computations are also 
chosen to be equal for each type component, and 

$1 = 3^, for i = 1, 2, . . . k. 

For purposes of comparison, it is assumed that each 
type of component is placed on test for an equal length of 
time, t^ = t^, and each type of component experiences the 
same number of failures during its total test time, 
fj^ = f^. This assumption is modified slightly to be able 
to examine the system that has only one component failure. 

A system exhibiting two component failures is also 
examined. 

The data (values of the parameters) have been chosen 
by the author. They have been intentionally chosen to be 
computationally simple yet to still exhibit the characteris- 
tics of a realistic system. 
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II. COMPUTATIONAL METHODS 



A. CLASSIC 

Reference 3 is used exclusively for the computations in 

this method. The computational formulas given in this 

reference are modified slightly to accommodate the basic 

assumptions of identical test time, identical number of 

failures for each component, etc. 

X. is a Maximum Likelihood Estimate for the failure 
1 

th 

rate of the i^ component. It will be equal to the total 
number of failures divided by the total test time, 
f . 

X = ^ 
i ti • 

X^, the system failure rate, will be equal to the sum 

of the individual component failure rates. X = 2 X. . 

Each component is assumed to fail independently and k w: .1 
equal the number of components that are placed together 
in a series. 



X = 
u 



2X^ + K’^C + (4Lk»^C + 



X^ will be an upper confidence level on the system 
failure rate. 



k X . 
C = ^ 
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The ’.’•alues for K’ are found in Table I. The computational 
procedure to determine the Beta values for a 90 per cent 
level of significance are also discussed. Eighty per cent 
level of significance values are found in Ref. 3. 

A formula for the upper confidence level for tke system 
failure rate is also provided when no failures have occurred 
during the total test time that each component has been 
allotted. That formula follows: 



A 



u 



K 



,2 



n 



k 

E 

i = l 




where n equals the number of component terms in the 
summation. 

Substitution of the upper limits on the failure rates 
obtained will generate corresponding lower confide'nce 
limits on reliability, thus 



= exp‘^u . 

Beta can be calculated from the following formula: 

X - X = Beta(K) (X 

X^ is obtained from Chi-Squared tables in Ref. 1 at 
the desired confidence level with 2(X+1) degrees of free- 
dom where X is the number of failures and K is the 
percentage point of the normal distribution point in the 
same confidence level. 
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TABLE I 



BETA VALUES FOR 90 PER CENT CONFIDENCE 



X 

No . of 
Failures 


^ .90,2(X+1) 


.9o'i2(X+l)^ 


Beta 

Value 


K‘ 

CBetaCl. 282)) 


0 


4.906 


2.3025 


1.18362 


1.5174 


1 


7.779 


3.889 


1.14271 


1.46496 


2 


10.645 


5.3225 


1.12336 


1.44015 


3 


13.362 


6.681 


1.1109 


1.42417 


4 


15.987 


7.994 


1.10188 


1.41261 


5 


18.549 


9.275 


1.09494 


1.403713 


10 


30.813 


15.406 


1.0743 


1.37725 


20 


54.090 


27.045 


1.05669 


1.35468 


30 


I N T 


E R P 0 L.A TED 




1.34639 


36.5 


91.1 


45.55 


1.04596 


1.34092 



The "Beta Value" listed in this table is germane to 
Ref. 3 and should not be confused with the 3 that is the 
optimistic scale parameter used in the Bayesian simulation. 

B. BAYESIAN 

The Bayesian technique will assume the apriori density 
to be Gamma ;a^ , 3^) . and 3^^ are the shape and scale 

parameters. The "prior" can be made "optimistic" if the 
scale parameter is chosen large in relation to the shape 
parameter . 
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The a posteriori density then becomes (X . ;a . +f - , 3 • +t . ) • 

t' 1 '•I’li’ii-' 

£. will equal the number of failures of the i^^ type com- 

th 

ponent and t^ will equal the total test time for the i*^ 
type component. 

k 

The distribution of X = Z X. will then be determined 

^ i=l ^ 

by computer simulation. Random variates of each X^ are 
generated from the Gamma distribution with parameters 
a. + f., 3- + t.. A random variate will be generated for 
each X^, i = 1, . . . k, the number of components in series. 
The X. , will then be added to determine the series system 
failure rate. 

The process for generating a system failure rate is 
then repeated 1000 times to yield ^s2’ ’ * ' ^slOOO' 

The 1000 random values of X are then ordered to yield 

s 

^s(l)’ ^s(2)’ • • • ^s(lOOO)* 

^su(y) "estimate" of the (1 - percentile 

point of the distribution of X^. 

The Bayesian 100(1 - y)^^ upper confidence level for 

Xs ^sl000(l - y)‘ 

At the time of this writing a subroutine to generate 
random variates from the Gamma distribution does not exist 
in the Naval Postgraduate School computer library. Reference 
2 was used to generate random variates from the Gamma 
distribution with an integer shape parameter. 
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Th-r . gejie rat ion of random variates WEth a shape parameter 
that is not an integer poses a muds, more complicated prob- 
lem. Reference 4 has been written to handle this situation 
for shape parameters between Q.05 and 1.0. Although 
LT Robinson's approach is relatively untried, it sho\irs a 
great deal of promise and may possibly be incorporated 
into the Naval Postgraduate School computer library in the 
future. It has been used in this writing for comparative 
purposes for a shape ‘parameter less than 1.0. 

C. SEMI-CLASSIC 

Similar to the Classic technique, Ref. 3 is used to 
compute the upper confidence levels for the system failure 
rate for this method. The computations will be modified 
somewhat by adding the identical "optimistic" priors used 
in the Bayesian technique. 

The number of failures per component will be added to 

and the total t^st time per component will be added to 

t h 

A maximum Likelihood Estimate for the i*^ component's 
failure is then: 



L ^ 



f . + a 



t. * e. 



The system failure rate for this Semi-Classic method is 
the same as for the Classic method. 






X 



s 
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The value of C will also change slightly: 



X. 

1 




When no failures have occurred, the upper confidence 
on the failure rate now is: 



2 

^ _ K'^ 1 

u n t. + 3. 

1 1 

The computation for the upper confidence level with 
remains the same. 

This method is somewhat "Bayesian" in the sense 
it is utilizing a prior but it is "Classic" because 
nature of the computational method; 



limit 



failure 

that 
of the 
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III. PARAMETERS OF VARIATION 



Although the range through which some of the parameters 
are varied could be more extensive, it does provide the 
reader with an adequate base for comparing the methods of 
computation. 



1 . 0 , 0. 5 

3. 50, 000 

t^ 30, 50, 100 mission units 

f^ 1.0, 0.0 The failures are also modified so 

that f^ = 1, f 2 = f^ = 0 and f^ = 1, f^ = 1, 

fj = f , fj^ = 0. This will provide the 

reader with the opportunity to compare a 
system where only one component or only two 
components fail. 



k 1, 2, 5, 10, 30 components 

Y 0.10 
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IV. COMPARISONS OF SYSTEM FAILURE RATE 



UPPER CONFIDENCE LEVELS 

The following 12 tables provide the reader with upper 
confidence levels computed from the same data by the three 
different methods: Classic, Semi-Classic, and Bayesian. 

A program to generate random variates from a Gamma 
distribution with a shape parameter that was greater than 
1.0 but not an integer was not available; therefore UCL's 
for the Bayesian simulation could not be computed for the 
case when the shape parameter was less than 1.0 and a 
failure existed for the i component. 

All calculations were conducted on the IBM 360 computer 
with the exception of the Classic method with zero failures, 
whose calculations were computed on a desk calculator. 

For each three-line block of numbers the reader may 
compare the three methods’ upper confidence levels com- 
puted from the same arguments. To see how a change in the 
scale parameter, the number of times a component fails, 
or the number of components in series is reflected in the 
upper confidence level, the reader merely moves to the right 
or down in the table. If the reader wants to see how a 
change in the shape parameter or a change in the amount of 
time a component is on test is reflected, he must go to 
another table. 
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SYSTEM FAILURE RATE UPPER CONFIDENCE LEVELS ALFA, 
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V, CONCLUSIONS 



A "crossover" point in this discussion will be defined 
as the point or general area at which an upper confidence 
level that was previously lower than one with which it was 
being compared becomes higher. 

The Bayesian method will be compared with the Classic 
and the Semi-Classic will be compared with Classic. The 
Bayesian and the Semi-Classic appear to behave similarly 
and remarks concerning this similarity will be mentioned. . 

A. UNMODIFIED COMPARISONS, ALFA =1.0 

The first three tables deal with a shape parameter of 
1.0 and a system of components that either all of the 
components experience a failure or none of the components 
experience a failure. 

In the cases where none of the components fail, a 
crossover point is exhibited in every case when five or 
more components are in series, regardless of how optimistic 
the scale parameter becomes or how long the test time is 
extended. When every component fails, the Bayesian and Semi- 
Classic systems do not crossover. In both cases the values 
of the upper confidence level for the Bayesian and the 
Semi-Classic methods become closer together as the number 
of components increases. 
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B. UNMODIFIED COMPUTATIONS, ALFA =0.5 



The upper confidence levels for the second three tables 
(less the Classic) are computed with a more optimistic 
prior shape parameter (alfa = 0.5) so naturally the UCL’s 
for these two methods are lower than those shown in the 
first three tables. 

Crossover points are much the same as in the case when 
alfa = 1.0. The Bayesian and Semi-Classic methods cross- 
over when none of the components fail; only the crossover 
point occurs when about 10 components are placed in series. 
The Semi-Classic does not crossover when all components 
fail (the Bayesian is not examined due to the non-existance 
of a Gamma random variate generator with a non- integer 
shape parameter) . 

C. MODIFIED COMPUTATIONS, ALFA =1.0 

The third three tables provide a much more 
interesting comparison. Here the situation where only 
one (or two) component (s) in the system fail. The Semi- 
Classic and the Bayesian techniques will crossover the 
Classic system in every case regardless of how optimistic 
the "priors” are or how long the components are tested. 

When the priors are more optimistic and/or the test time is 
larger, then it takes a larger number of components in 
series for the crossover to occur. This point may occur 
with as few as two components or as many as between 10 and 
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30 components. The Bayesian UCL's are strictly lower than 
the Semi-Classic UCL's. 

D. MODIFIED COMPUTATIONS, ALFA =0.5 

In the last three tables only the Semi-Classic and the 
Classic methods can be compared for the Bayesian cannot be 
simulated for this case (non-integer). 

The Semi-Classic technique is using a more optimistic 
prior so its UCL's are lower than before but again in every 
case it will 'crossover with the Classic method. 

Generally one could conclude that as long as the 
number of components that are in series in a system are 
few, the Bayesian approach may yield a lower value of an 
upper confidence level for the system failure rate, pro- 
vided the priors can be chosen optimistically. If the 
system is complex enough that more than "a few" components 
must be placed in series then the Classic system appears to 
yield the lowest values of upper confidence levels. The 
Semi-Classic technique would only be useful if the priors 
were optimistic, the number of components in series were 
few, and the optimistic priors were non- integers . 



/ 
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COMPUTER PROGRAM FOR THE CLASSIC METHOD 



THE VARIABLE DEFINITIONS FOR THIS PROGRAM FOLLOW: 

TIME THE NUMBER OF MISSION UNITS EACH COMPON- 

ENT IS UNDER TEST 

FAIL THE NUMBER OF TIMES THAT EACH COMPONENT 

WILL FAIL DURING THE TOTAL TEST TIME 
XK THE NUMBER OF COMPONENTS THAT ARE IN THE 

SERIES SYSTEM 

XLAMI INDIVIDUAL COMPONENT FAILURE RATE 

XLAMS SYSTEM FAILURE RATE 

XKPRI A VALUE TAKEN FROM TABLE I (BETA VALUES) 



XK=1.0 

XKPRI=1 .469496 
7 TIME=30.0 

10 C=0.0 

FAIL=1.0 

XLAMUP=0.0 

XLAM I=FA IL/TIME 

XLAMS=XK*XLAMI 

C = ( XK^ { XLAMI /T I ME) ) /XLAMS 

XLAMUP= ( 2.0 -'XL AMS+( XKPRI )=^C + SORT( ( (4. XLAMS) =?'( XKPR 
1 I«^2 )>=C )-MXKPRI^^4)4'.C-*2) )/2.0 
WRITE! 6, lODFAIL ,XK,TI MEtXLAMUP 

101 FORMAT!' • , 5X , ’ F A I L= ' , F7 .3 » 5X , • XK= ’ , F 7 . 3 ♦ 5 X , • TI ME = • , F 7 
2.3,5X, • XLAMUP = ',F10.6,/) 

TIME=T1ME-^20.0 

102 IF(T1M^.LE.51.0) GO TO 10 
TIME=TIMF+30.0 

103 IFITIME.LF.IOI .0) GO TO 10 
XK=XK-H.O 

XKPR1=1. 44015 

104 IFIXK.LE.2.3) GO TO 7 
XK = XK-f2.0 

XKPRI=1 .403713 

105 IF(XK.LT.5.5) GO TC 7 
XK = XK-i-2.0 

XKPRI=1 .37725 

106 IF!XK.LT.10.5) GO TC 7 
XK=XK-H5.0 

XKPRI=1 .34092 

107 IF! XK.LT.30.5) GO TC 7 
STOP 

END 
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COMPUTER PROGRAM FOR THE SEMI -CLASSIC METHOD 



THIS SYSTEM IS LABELED "SE M I-C L A S S I C '• BECAUSE IT 
UTILIZES THE "PRIORS" THAT WERE INPUTS TO THE BAYES- 
IAN TECHNIQUE BUT THE COMPUTATIONS ARE PERFORMED THE 
THE SAME WAY AS IN THE "CLASSIC" TECHNIQUE. 

THE VARIABLE DEFINITIONS FOR THIS PROGRAM FOLLOW: 

XK 

FAIL 
TIME 

BETA 
ALFA 
XLAMS 
XLAMI 
XLAMUP 



ALFA=0.5 
XK=1. 0 
XKPR1=1 .469496 
FA1L=0.0 
BETA=50.0 
TIME=30.0 
0 C=0.0 

XLAMUP=0.0 

XLAMI= ( ALFA+FAIL )/ (TIME+BETA) 

XLAMS=XK*XLAMI 

C={ XK«( XLAMI/( TIME+BETA) ) ) /XLAMS 

XLAMUP = ( (2 .0-^XLAMS ) + ( ( XKPR I’^=«2 )*C ) + SORT( { t 4.0*XLAMS) 
1*( XKPRI-«=*2)^'C) + (XKPRI*-=^4)*C>^*2) )/2.0 
WRITE! 6, 10 DAL FA, FAIL, XK , T I ME , B E TA , XL AMUR 

101 FORMAT!' • ,2X , ' ALFA=' , F4. 1 ,2X , • FA IL= ' , F4. 1 , 2X, ' XK= • 
6F5.1,2X,'TIME=' ,F6.1,2X,'BETA=* , F6 . 1 , 2X , • X L AMUP= ' , 
7F10.6,// ) 

TIME=TIME+20.0 

102 IF!TIME.LE. 51.0) GO TO 10 
TIME=TIME+30.0 

103 IF!TIME.LE . 101 . 0) GO TO 10 
BETA=3ETA+50.0 

104 IF !BETA. LE. 101 .0 ) GO TO 9 
FAIL=FAIL+1.0 
IF!FAIL.LE.1.1 ) GO TO 8 
XK=XK+1.0 

XKPRI=1. 44015 

IF!XK.L5.2.1 ) GO TO 7 

XK=XK+2. 0 

XKPRI=1. 403713 

IF!XK.LE.5.1 ) GO TO 7 

XK=XK+2.0 

XKPRI = 1. 37725 

IF!XK.LE.10. 1) GO TO 7 

XK=XK+15.0 

XKPRI=1 .34092 

IF !XK. LE.30.5) GO TO 7 

ALFA=ALFA+0. 5 

IF!ALFA.LE.1.1 ) GO TO 6 

STOP 

END 



NUMBER CF COMPONENTS IN SERIES 
NUMBER CF FAILURES PER COMPONENT 
NUMBER OF MISSION UNITS EACH COMPONENT 
IS UNDER TEST 

OPTIMISTIC SCALE PARAMETER 
OPTIMISTIC SHAPE PARAMETER 
SYSTEM FAILURE RATE 

ESTIMATE OF THE COMPONENT FAILURE RATE 
UPPER CONFIDENCE LEVEL ON SYSTEM 
FAILURE RATE 
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COMPUTER PROGRAM FOR THE MODIFIED SEMI-CLASSIC METHOD 



THIS PROGRAM IS A "MCDIFIEO S EMI -CLASS I C" IN THAT 
THE FAILURES ARE MODIFIED SO THAT ONLY ONE COMPONENT 
OR ONLY TWO COMPONENTS FAIL DURING THE COMPLETE SYS- 
TEM TEST. ALL OTHER VARIABLES ARE THE SAME AS FOR 
THE "SEMI-CLASSIC" PROGRAM. 



ONLY ONE set OF VARIABLES ARE TESTED WITH THIS 
PROGRAM. THE VARIABLES T I M E , B ET A , AL FA , AND FAIL 
WOULD HAVE TO BE VARIED TO GENERATE A FULL TABLE OF 
VALUES. 



XKPRI=1 . A4015 
XK=2.0 
TIME=10D.0 
XLAMS=2. 0/100.0 
BETA=0.0 
10 ALFA=0.0 

C= 1.0/ TIME 
XLAMUP=0.0 

XLAMUP= ( 2. 0«XL AM S+( XKPR 1**2) '-"0 + SORT ( ( ( 4 . 0* XLAMS ) * ( XKPR 
1 I**2)*C)4(XKPRI**4)*C**2) )/2.0 
WRITE! 6, 101) XK .BETA ,TI ME , XL AMUR 
101 FORMAT!' •,2X,'XK = • , F7 . 3 , 2X , ' BE TA = ',F7.3,2X, 

2'TIME = • ,F7.3,2X, 'XLAMUP = '.FIO.B,/) 

XK=XK+3.0 

XLAMS = 2 .0/100. 0 

IFIXK.LE.5.5) GO TO 10 

XK=XK+2.0 

XLAMS=2. 0/100.0 

IFIXK. LE. 10. 5) GC TO 10 

XK=XK+15.0 

XLAMS=2 .0/100. 0 

IF(XK.LE.31.0) GO TO 10 

STOP 

END 
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COMPUTER PROGRAM FOR THE BAYESIAN SIMULATION 
ALFA = 1.0 



THE DEFINITIONS CF THE VARIABLES USED IN T^HIS 
PROGRAM FOLLO»i: 



BETA 

ALFA 

FAIL 



K 



XLAMS 



TIME 



XLAM 



THE NUMBER OF MISSION UNITS THAT EACH 

COMPONENT IS UNDER TEST 

THE OPTIMISTIC SCALE PARAMETER 

THE OPTIMISTIC SHAPE PARAMETER 

THE NUMBER OF TIMES THAT EACH COMPONENT 

WILL FAIL DURING THE TEST TIME 

THE NUMBER OF COMPONENTS THAT ARE IN THE 

S ER I ES S YST EM 

THE FAILURE'rATE FOR AN INDIVIDUAL 
COMPONENT 

THE SYSTEM FAILURE RATE 



THIS SYSTEM IS CALLED 
DIMENSION XLAMS (1000) .KEY(IOOO) 

TIME=30. 0 
^ BETA=50.0 
5 FAIL=0.0 
7 IND = 1 

600 DO 610 10=1. 1000 
XLAMS(IO) = 0.0 
610 CONTINUE 
6 ALFA=1.0 

B=TIME+BETA 
KA=ALFA + FAIL 

GO TO (601,602.603.604.606.607) . IND 
6 01 K = 1 

GO TO 608 

602 K = 2 

GO TO 608 

603 K = 5 

GO TO 608 

604 K = 10 

GO TO- 608 
606 K = 30 

608 IX=999 

9 DO 50 J=1.1000 

DO 609 ILBDS = l.K 
TR=1.0 

10 DO 20 1=1. KA 

CALL RANDU( IX. lY.YFL ) 

IX = I Y 

20 TR=TR^YFL 
30 XLAM=-ALOG(TR)/B 

XLAMS(J) = XLAMS(J) + XLAM 

609 CONTINUE 
KEY( J)=J 

50 CONTINUE 

CALL SHSORT( XL AMS. KEY. 1000) 

WRITE(6.60 )K.TIME, B ET A . FA I L . XL AMS ( 900 ) 

60 FORMAT!' '.3X.'C0MP. = • . I 4 . 3X . • T I ME = '.F7.2.3X. 
I'BETA = '.F7.2.3X. 'FAIL = • . F6 . 2 . 3X. • LA MS = '. 
2F10.6. //) 

IND = IND + 1 
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GO TO 600 
6C7 CONTINUE 

FAIL=FAIL+1-0 

IF(FAIL .LE.l .5 ) GO TO 7 

3ETA=BETA+5C.O 

IF(BETA.LE. 100.1) GO TO 5 

TIME=TIME+2C.O 

IF( TIME. LE . 50. 1) GO TO 4 

TIME=TlME+30.0 

IF(TIME.LE. 100.1) GO TO 4 

STOP 

END 
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COMPUTER PROGRAM FOR THE 
MODIFIED BAYESIAN SIMULATION 
ALFA = 1.0 



THIS BAYESIAN TECHNIQUE IS ’'MODIFIED" IN MUCH THE 
SAME WAY THAT THE SEMI-CLASSIC TECHNIQUE WAS IN THAT 
THE NUMBER OF FAILURES ARE "MODIFIED" SO THAT ONLY 
ONE AND ONLY TWO COMPONENTS WILL FAIL DURING THE SYS- 
TEM TEST TIME 



THE DEFINITIONS OF THE VARIABLES FOR THIS PROGRAM 
FOLLOW: 

TIME THE NUMBER OF MISSION UNITS THAT EACH 

COMPONENT IS UNDER TEST 
BETA THE OPTIMISTIC SCALE PARAMETER 

K THE NUMBER OF COMPONENTS THAT ARE IN 

THE SERIFS SYSTEM 

JJ THE NUMBER OF COMPONENTS IN THE SYSTEM 

THAT FAIL 

XLAMS THE SYSTEM FAILURE RATE 



DIMENSION XLAMS (1000) , KEY (1000) 

1 DO 700 JJ=1,2 
TIME=30.0 
A BFTA=50.0 
5 TND=1 

KA IS EQUAL TO ALFA PLUS FAIL» WHICH IS 2 
THE INDEX "JJ" IS EQUAL TO THE NUMBER OF COMPONENTS 
THAT FAIL IN THE SERIES SYSTEM 

600 DO 610 10=1 ,1000 
XLAMS( I0)=0.0 

610 CONTINUE 

B=BFTA+TIME 

GO TO ( 601, 602, 6C3t 60At606,607) ♦ IND 

601 K=1 
GO TO 608 

602 K=2 
GO TO 608 

603 K=5 
GO TO 608 

604 K=10 
. GO TO 608 

606 K=30 
608 IX=999 

DO 80 J=l,1000 
DO 24 JKJ=1,JJ 
TR = 1 .0 

10 DO 20 1=1,2 

CALL RANDU( I X, I Y,YFL) 

IX = IY 

20 TR=TR*YFL 

XLAMX=-ALOG(TR) /B 
XLAMS (J )=XL AMS (J )+XL AM X 
24 CONTINUE 

IF(K-JJ .LE.O) GO TO 79 
L=K-JJ 

DO 30 ILBDS=1,L 
TR=1.0 

C KA IS NOW EQUAL TO ONE FOR THE REST OF THE SYSTEM 

DO 25 I I =1 , 1 
CALL RANDU( I X, lY, YFL ) 
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IX = IY 

25 7R=TR*YFL 

XLAM=-ALOG(TR)/B 
XLAMS( J)=XLAMS(J )+XLAM 

30 CONTINUE 

79 KEY(J)=J 

80 CONTINUE 

CALL SHSORTI XLAMS, KEY, 1000) 

WR IT E( 6, 90 )K, TIME, BETA, XL AMS (900) 

90 FORMATC ',3X,'CCVP. = • , 1 4 , 3X , • T I ME = ',F7.2,3X, 
*'BETA = • ,F7.2,3X, ' XLAMS = *,F10.6,//) 

IN0=IN0+1 
GO TO 600 
607 CONTINUE 

BETA=BET A+50.0 
IF(3ETA.LE. 100.1) GO TO 5 
TlME=TIME+20.0 
IF(TIME.LE.50.1 ) GO TO 4 
T1ME=TIME+30.0 
IF(TIME.LE.100.1 ) GO TO 4 
7C0 CONTINUE 
STOP 
END 
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COMPUTER PROGRAM FOR THE BAYESIAN SIMULATION 
ALFA =0.5 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 



THIS PROGRAM WILL GENERATE VARIATES FROM THE GAMMA 
0ISTRI8UTIGN WHEN ThE SHAPE PARAMETER (ALFA) IS LESS 
THAN l.D 



THF VARIABLE DEFINITIONS FOR THE MAIN PROGRAM FOLLOW 



K 

XBETA 

XTIME 



XLAMS 

XLAMUP 



XFAI L 



NCMP 



TFF OPTIMISTIC SHAPE PARAMETER 
THE OPTIMISTIC SCALE PARAMETER 
THE NUMBER OF MISSION UNITS THAT EACH 
COMPONENT IS UNDER TEST 

THE NUMBER OF TIMES THAT EACH COMPONENT 

FAILS CURING THE TEST TIME 

THE NUMBER OF COMPONENTS IN THE SERIES 

SYSTEM 

THE SYSTEM FAILURE RATE 

THE UPPER CONFICENCE LEVEL ON THE 

FAILURE RATE 



DIMENSION XLAMSI 1000) ,KEY( 1000) 

REAL^4 K 
K = 3. 5 

XTIME=30.0 
30 XBETA=50.0 
40 IN0=1 

600 CO 610 10=1,1000 
XLAMS( I0)=0.0 

610 CONTINUE 

BETA=1 .0/(XBETA+XTIME) 

GO T0( 601,602,603, 604, 605, 606) , IND 

601 NCMP=1 

GO TO 607 

602 NCMP =2 

GO TO 607 

603 NCMP=5 

GO TO 607 . 

604 NCMP=10 
GO TO 607 

605 NCMP=30 
607 IX=999 

CALL GMINITIK, BETA) 

10 DO 50 J=l, 1000 

CO 609 ILBCS=1 , NCMP 
CALL GAMAI IX,Z) 

C Z WILL = AN UNORDEREC VALUE OF LAMBDAI 

XLAMSI J )=XLAMS ( J )+Z 
609 CONTINUE 
KEY( J)=J 
50 CONTINUE 

C NOW ORDER THE SYSTEM FAILURE OATES 

CALL SHSORTI XLAMS, KEY, 1000) 

WRITE (6 ,60 ) NCMP, XBETA, XTIME, XL AMS (900) 

60 FORMATI' ',2X,'NCMP = ' , 1 4 , 2X , • B ET A =',F6.2,2X, 
>^'TIME =' ,F6.2, 2X, ' XLAMS =',F10.6,/) 

IND=IND+1 
GO TO 600 
6 06 CONTINUE 

XBETA=X3ETA+50.0 
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IF(XBETA.LE. 100.1) GO TO 40 
XTIME=XTIME+20 .0 
1F(XTIME.LE.50.1) GO TO 30 
XTIME=XTIME + 30 .0 
IF(XTIME.LE.100.1 ) GO TO 30 
END 
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SUBROUTINE GMIN I T ( K , B ET A ) 
ADDITIONAL ENTRY POINTS: 
RSULT 

GAMA( IX, Z) 



PURPOSE: 

GENERATION OF GAMMA RANDOM DEVIATES WITH SHAPE 
parameter less than ONE. 



METHOD: 

A MODIFICATION OF MARSAGLIA'S BCX-W EDGE-T A I L METHOD 
FOR NORMAL DEVIATES IS USED. THE PDF IS DECOMPOSED 
INTO A HEAD REGION, A NUMBER (DEPENDENT ON K) OF 
RECTANGLES AND WEDGES AND A TAIL REGION. THE GMINIT 
SECTION OF THE SUBROUTINE ALSO SETS UP A BINARY 
SEARCH TREE TO BE USED FOR EFFICIENT SELECTION OF THE 
PROPER REGION DURING THE ACTUAL GENERATING PROCESS, 
WHICH IS HANDLED BY THE GAMA SECTION 

DESCRIPTION OF PARAMETERS 

K GAMMA DISTRIBUTION SHAPE PARAMETER 

(MUST BE .GE. 0.05 AND .LE. 1.0) 

BETA GAMMA DISTRIBUTION SCALE PARAMETER 

IX SEED FOR RANDCM NUMBER GENERATOR 

Z RETURNED GAMMA DEVIATE 

THE PDF OF THE GAMMA f^UNCTION IS GIVEN BY 

F(X) = (1/BETA)=^*K ^ X**(K-1) * E XP (- X/B E TA ) /GAMMA ( K ) 

THE FOLLOWING SUBROUTINES ARE USED: 

RANDOM! IX) RETURNS A UNIFORM (0,1) DEVIATE 
INVGAM(K,X) COMPUTES THE INVERSE GAMMA CDF 
IGAM(K,X) COMPUTES THE INCOMPLETE GAMMA 
FUNCT ICN(GAMMA CDF) 



NOTE: 

UNDERFLOW IS POSSIBLE WHEN K IS LESS THAN. 18 AND 

BECOMES MORE LIKELY AS K DECREASES. WHEN K IS 0.5 
THE PROBABILITY 0^^ UNDERFLOW IS ABOUT. 0C0129 FOR 
ANY GENERATED DEVIATE. 



SUBROUT INE 

REAL*4 K,I 

INTEGER*4 

L0GICAL=^1 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

EOUIVALENC 



GMI NIT( K 
NVGAM. IGA 
FIRST, TAB 
USED 

P( 100) ,X( 
TABLE ( 202 
TEST (202 ) 
RAND( 2) 

E (U.RAND 



.BETA) 

M 

LE, BOTTOM, END 

101) ,H( 100) ,0(100) ,R 
), PROB( 202) , NEXT! 202 
,LIST(202 ) ,USED(202) 

(1 ) ), (V,RAND(2) ) 



(100) ,B(100 
) ,LAST(202) 



) 



THIS FIRST SECTION 
BE USED BY GAMA 
ARE USED BY GAMA: 

PO PROBABILITY 

PN PROBABILITY 

P(I) PROBABILITY 

X(T) LEET-HAND B 

H( I ) WIDTH OF I- 

0(1) PROBABILITY 



INITIALIZES CONSTANTS AND TABLE 
WHEN IT IS CALLED. THE FOLLOWI 



FOR "HEAD" REGION 
FOR "TAIL" REGION 
FCR I-TH RECTANGLE 
CUNDARY OF I-TH RECTANGLE 
TH RECTANGLE 
CF I-TH WEDGE 



S TO 
NG 
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R(I) REJECTION TEST RATIO FOR I-TH WEDGE 
B(I) Y INTERCEPT FOR I-TH WEDGE 
ALPHA SHAPE PARAMETER - 1. 

PROS ORDERED VECTOR CF PROBABILITIES 
TABLE VECTOR OF W E DGE /R EC TANGL E NUMBERS CORRESPOND- 
ING TO PROBABILITIESS IN PROB 
FIRST STARTING POINT FOR BINARY SEARCH 
NEXT, LINKS FOR BINARY SEARCH 
LAST 

J1 POSITION IN PROB OF PO 

HI TO CONSTANTS FOR APPROXIMATION TO INVERSE GAMMA 
HA CDF FOR SMALL VALUES OF Z 



CHECK FOR K IN RANGE 

IF( (K.GE.0.05) .AND. ( K . LE . 1 . 0 ) ) GO TO 5 
WRITEI6 ,4) K 

4 FORMAT! //• OGMI NIT CALLED WITH K=',1PE16.6, 

OUT OF RANGE'/) 

RETURN 

GET UPPER BOUND ON NUMBER OF RECTANGLES 

5 N=20.+6.6/K 
IFIN.GT . 100)N=100 
M = 2'^N+2 

MM=M-1 

ALPHA=K-1. 

GK=GAMMA (K ) 

P0=5.E-5/(K=^K) 

HFAC=2. 

SET UP RECTANGLE BOUNDS 

X(1)=INVGAM(K, PO) 

P0=IGAM(K,X(1) ) 

H( 1 )=.25^X { 1 ) 

CO 10 I=2,N 

X( I )=X(I-1)+H{I-1) 

H( I ) = H( 1-1 )VHFAC 

pm=o. 

0( I ) = 0. 

10 CONTINUE 

X(N+1) =X(N)+H(N) 

ZERO PROBABILITY VECTORS AND LINKS 
DO 15 I =1,M 

NEXT! I )=0 
LAST! I )=0 
PROB! I ) =0. 

LIST!! )=0 
USED! I )=. FALSE . 

15 CONTINUE 

FIND COEFFICIENTS FOR NEWTON-R APHSON APPROXIMATION TO 
SLOPE OF PDF 

B1=-2.«ALPHA 
B2=ALPHA-« ! ALPHA-1 . ) 

A1=B1 + 1 . 

A2=ALPHA^! ALPHA-2. ) 

C=l. -ALPHA 

FIND RECTANGLE PROBABILITIES AND WEDGE VALUES 
PL=P0 

FL = EXP!ALPHA='-AL0G(X!1) )-X! 1) )/GK 
DO 40 1=1, N 

FU = EXP! ALPHA=^-ALOG!X! I + l ) )-X! I + l ) ) /GK 
PU = IGAM!K ,X ! I + l ) ) 

P ! I )=H! I )^FU 
0! I ) = PU-PL-P! I ) 
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NEWTON-RAPHSON ITERATION TO FIND POINT WHERE 
SLOPE OF PDF IS PARALLEL TO CHORD 

W=X(I ) 

S={ FU-FL) /H( I ) 

SC=S*GK 
DO 20 J=l,8 

Y=W*( (W+Al)*W+A2+SC*EXP(C*AL0G(W)+W))/( (W+Bl 
IF(A8S(Y-W) .LT.1.E-4=^Y)G0 TO 30 
W = Y 

20 CONTINUE 

FIND VALUES FCR REJECTION METHOD 

30 B( I )=EXP( AL PH ANALOG (Y)-YJ /GK+S*( X( I )-Y) 

R( n = (B(I )-FU) /(FL-FU) 

TEST TO SEE IF ENOUGH RECTANGLES HAVE BEEN TAKEN 
IF( PU.GT. 0.999)GO TO 45 

RESET PROBABILITY AND FUNCTION VALUES FOR NEXT 
RECTANGLE 

PL = PU 
FL = FU 

40 CONTINUE 

FIND LOWER END OF NON-ZERO PROBABILITY VECTOR 
45 L0W=2*( N-I )+l 

FIND TAIL PROBABILITY 

SUM1=P0 

DO 50 J=1.N 

SUM1=SUM1+P( J)+0( J) 

50 CONTINUE 
PN=1.-SUM1 

GET VALUES TO COMPUTE LOW END APPROXIMATION TO INVERSE 
GAMMA CDF 

H1=K*GK«P0 

H2=l./K 

H3=(K+1 . )Y5.0E-7 
H4=4./(K+1. ) 

GENERATE PROBABILITY VECTOR 

DO 80 1=1, N 

PR0B( I )=P ( I ) 

TABLE ( n = I 
PR08( 1+N) =0(1 ) 

TABLE! I+N )=-I 
80 CONTINUE 

PROB(M-1)=PO 
TABLE( M-1) =0 
PR0B(M)=PN 
TABLE(M) =0 

SORT PROBABILITY VECTOR 

DO 110 1=1, MM 
ICH=0 
L = M-I 

DO 100 J=1,L 

IF(PR0B( J )-PROB( J + 1) ) 100, 100,90 
90 TEMP=PROB(J) 

PROB ( J) =PR0B ( J+1) 

PROB( J + 1 )=TFMP 
ITEMP=TAELE(J) 

TABLE! J) =TA8LE( J+1) 

TABLE! J + 1 ) = ITEMP 
ICH=1 
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100 

110 



120 



130 

140 

150 
1 51 

152 

153 

154 

155 



156 

157 

1571 

158 
1585 

159 

1591 

160 

161 

162 

163 



CONTINUE 

1F( ICH) 120, 120,110 
CONTINUE 
PROS (M) = 1. 

CONVERT PROS TO COPULATIVE PROBABILITIES AND 
FIND FIRST AND J1 

J1 = 0 
FIRST=0 
L=LOW+l 
DO 130 I=L,M 

IF( (TABLE ( I ) .EO.O) .AND. ( PROB { I ) . EO. PO) ) J1 = I 
PR03( 1 )=PR0B( I )+PR0B(I-l) 

IF( (PRCB{I).GE.D.5) .AND. ( F I RST . EQ.O ) ) F I RST= I 
CONTINUE 

IF(FIRST.EO.M) FI RST=MM 

IF(Jl) 14C, lAO, 150 
J1=L0W 

NOW DETERMINE THE VECTORS NEXT AND LAST FOR BINARY 
SEARCH 

B0TT0M=1 
END = 1 
PR=.25 

LIST ( 1 ) =FI RST 
TEST( 1 ) =.5 
USED(FIRST)=.TRUE. 

DO 159 1=1, BOTTOM 
LI =LI ST(I ) 

I F(USEC(L I+l) )G0 TO 155 
PRN=TEST( I )+PR 
DO 152 J = LO'a,MM 

IF(PROB( J) .GT.PRNIGO TO 153 
CONTINUE 

IF(.NOT. USED(J))GC TO 154 
J = J-1 

IF ( LI-J )153 ,155,155 
NEXT! LI )= J 
L IST( I )=N5XT(L I ) 

IF( (LI.EO.LOW) .OR. US ED { L I-l ) ) GO TO 159 
PRL=TEST( I )-PR • 

TEST! I ) = PRN 
DO 156 J=LOW,MP 

IF(PROB( J) .GT.PRDGO TO 157 
CONTINUE 

IE(.NOT. USED(J))GC TO 1571 
J=J+1 • 

IF(LI-J)158,158,157 
•LAST{ LI ) = J 
GO TO 153 5 
J = LI 

EN0=END+1 
LIST( END)=J 
TEST( END) =PRL 
CONTINUE 
BOTTCM=END 
PR = PR^O. 5 
1 = 1 

IF(LIST( I) )160,160,163 
B0TTCM=B0TTQM-1 
IF(BOTTOM) 165,165, 161 
DC 162 J=I , BOTTOM 

LI ST( J) =LI ST{ J+1) 

TFST( J)=TEST( J+1) 

CONTINUE 
GO TO 1591 

USEO(LIST( I ) )=.TRUE. 

1=1 + 1 

IF( I .LE.BOTTOMIGO TO 1591 
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EKD=BOTTCM 
GO TO 151 

SETUP FIRST CALL TO RANDOM 
165 CALL OVFLOW 

THROUGH WITH SETUP PHASE. QUIT. 
RETURN 



THIS SECTION PRINTS OUT THE VALUES GENERATED BY GMIMIT 

ENTRY RSULT 
WRITF(6,166)K, BETA 

166 FORMAT! /' IGENERATED VALUES FOR K=',1PE14.6» 

BETA=* ,E14.6) 

WRITE(6, 170) PO .PN 

170 FORMAT! • OHEAD PROB A BI L I T Y= ' , IP E 1 5 . 6 , 

TAIL PRG3A3ILITY = * , E15.6 ) 

WRITE! 6. 180) 

180 FORMAT! ' ORECTANGLE/WEDGE VA LUE S * //2X , ' 1 ' , 9 X , ' X ! I ) ' , 1 2X 
^ P! I ) • ,12X, '0 ! I ) ' , 12X , 'R! I ) • ,12X, ' B! I ) '/) 

SUM1=0. 

SUM2=0. 

CO 192 1=1 .N 
IF!P! I ) ) 193,193, 185 
185 SUM1=SUM1+P! I ) 

SUM2 = SUM2+0! I ) 

WRITE! 6, 190) I, X! I ) , P! I ) ,0(1 ) ,R! I ) ,B! I ) 

190 F0RMAT!1XI3,1P5E16.6) 

192 CONTINUE 

193 WRITE! 6, 194) SUMl ,SUM2 

194 FORMAT! 'OSUMS' , 1 5X , 1 P 2 E 16 . 6 ) 

WRITE! 6, 195) J1 ,H1 ,H2-,H3,H4 

195 FORMAT !/ 'OVALUES FOR HEAD/TAIL APPROXIMATION:*/' Jl=', 

H1 = ',E16.6,' F2 = ',E16.6,' H3 = ',E16.6,' H4', 

1E16.6) 

WRITE!6.196)FIRST 

196 FORMAT! / 'OSTARTl NG POINT FOR BINARY SEARCH', 14) 

WRITE ! 6, 197) ! I , PROB! I ) , TABLE! I ) ,NEXT! I ) , LAST! I ) ,I =1 ,M) 

197 FORMAT !/ '0PR03ABIL ITY VECTOR AND L INKS '// 14X, ' PROB ', 9X 
. 'TABLE ' ,4X,'NEXT',5X,'LAST'//!1XI5,1PE16.6,0P3I9)) 
RETURN 

THIS IS THE SECTION THAT ACTUALLY COMPUTES THE RANDOM 
VARI ATES 

ENTRY GAMA! IX, Z) 

GET TWO UNIFORM DEVIATES 
CALL RANDOM! IX, U, 2) 

CONDUCT BINARY SEARCH OF PROBABILITY VECTOR 
J=F1RST 

200 IF!U-PR0B!J ) )210,25C,230 
210 IF!LAST!J) ) 250,250,220 
220 J=LAST!J) 

GO TO 200 

233 IF!NEXT( J) ) 245 ,245 ,240 
240 J = NEXT ! J ) 

GO TO 200 
245 J=J+1 

LOCATED PROBABILITY DIVISION. GET TABLE VALUE. 

250 N=TABLE!J) 

IF!N)260,290,320 

THIS SECTION IS FOR T^E WEDGES 

260 N=-M 

CALL RANDOM! IX,U,1) 
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270 IF(U.LE.V)G0 TO 280 
TFMP=U 
U=V 

V=TEMP 

280 Z=X(N) +H(N )*U 

IF( V.LF.RC N) )G0 TO 330 

THIS STEP IS PERFORMED VERY RARELY 
W=U+EXP( ALPHA4:ALCG ( Z)-Z) /B(N) 

IF (V.LE.W) GO TO 330 
CALL RAMD0M( IX ,U,2 ) 

GO TO 270 

THIS SECTION IS FOR HEAD/TAIL PROBABILITIES 

290 IF( J.EO.Jl )G0 TO 300 
Z=INVGAM(K ,1 .-PN*V ) 

GO TO 330 
300 Z=(H1^V)==«=Fh2 

IF(Z.LT.H3)G0 TO 330 
l = l. + S0RT( l.-HA’^'Z) ) 

GO TO 330 

THIS SECTION IS FOR THE RECTANGLES 

320 Z=X( N) +H (N )^V 

SCALE GAMMA VARIATE AND EXIT 

330 Z=Z*BETA 
RETURN 
END 
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INVGAM SUBROUTINE 



FUNCTION INVGAM(K,Z) 

RFAL=^4 INVGAM, IGAM,K, EPS/1. E-08/ 

IF(Z.GT.O.O)GO TO 10 

INVGAM=0.0 

RETURN 

X=Z 

t=K-l. 

G=GAMMA ( K) 

DO 40 1=1,30 

IGAM(K,X)-Z)/ (EXP(T«ALOG(X)-X) ) 
IF(Y.LT.X)GO TO 30 
Y=Y=^.5 
GO TO 20 
X = X— Y 

IF( ABS ( Y)-EPS*X ) 50, 5C, 40 

CONT INUE 

INVGAM=X 

RETURN 

THIS IS THE END OF THE INVGAM SUBROUTINE 
END 
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IGAM SUBROUTINE 



FUNCTION IGAM(K.X) 

IMPLICIT REAL=!'8 (D) 

REAL*4 IGAM,K 
REAL^^B EPS/1. D-13/ 

IF( X.GT.O. »GO TO 5 
IGAM=0. 

RETURN 

5 IF(X.LE. 12.0)G0 TO 8 
IGAM=1.0 
RETURN 

8 DX=08LF(X) 

CK=DBLE( K) 

DTERM=DX**DK/DGAMMA (CK ) 

DSUM=DTERM/CK 
CO 20 1=1,30 

IF (DABS ( DTERM)-EPS=^DSUMJ A0,4Q,10 
10 DTERM=DTERM^(-OX )/DFLOAT( I ) 

CSUM=DSUM+DTERM/ (DK+CFLOAT ( I ) ) 

20 CONTINUE 
40 IGAM=SNGL( DSUM ) 

RETURN 

C THIS IS THE END OF THE IGAM SUBROUTINE 

END 
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